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The copper-oxygen layers of the parent compounds of the high-Tc superconductors are 
good physical realizations of the two-dimensional (2D) antiferromagnetic Heisenberg model 
. Recently, significant progress has been made in the theory of clean 2D antiferromagnets 
with short-range interactions [0,0. The nonlinear a-model (nlam) in 2 + 1 dimensions is 
believed to describe their long wavelength physics. Chakravarty et al. concluded that 
the ground state of this field-theory can be either ordered or disordered, depending on the 
coupling constant which depends on the interactions and the magnitude of the spin. There 
is convincing numerical evidence that the S = \ Heisenberg model with nearest-neighbor 
interactions on a square lattice is on the ordered side of the phase diagram, with a sublattice 
magnetization m ^ 0.3 0], which is close to the ordered moment observed in the cuprates 
i- 

The microscopic mechanism leading to the destruction of the antiferromagnetic long- 
range order upon doping the cuprates is still not understood. In particular, in spite of 
some experimental support it is not clear whether the nlam indeed describes the 

phase transition, as the dopants are expected to induce strong randomness in the spin- 
spin interactions before the charge carriers become mobile 0. It is therefore important to 
investigate how various types of quenched disorder affect the long-range order of the 2D 
Heisenberg model, and how the nlcrm description is altered by randomness. 

In this Letter we present quantum Monte Carlo results for the ^ = | random exchange 
model 

H =Y^ JijSi ■ Sj, (1) 

(id) 

where is a pair of nearest-neighbor sites on a 2D square lattice. The couplings Jij 

take two values, Jij = J(l ± A), at random, with a probability p for 1 + A and 1 — p for 
1 — A. We consider only the case A < 1, i.e. all couplings are antiferromagnetic and the 
system is non-frustrated. Although the hamiltonian does not exactly represent the kind 
of disorder present in the cuprates, where one would expect light doping to cause random 
frustrated interactions (corresponding to A > 1), we expect it to be of relevance in gaining 
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understanding of the effects of randomness in quantum antiferromagnets. In 2D an order- 
disorder transition can only take place at zero temperature, where the critical behavior is 
governed by quantum fluctuations. The opportunity to numerically study a disorder-driven 
quantum phase transition is a further motivation for this work. 

We study the behavior of the hamiltonian (|I|) in the {p, A) plane. For A — > the clean 2D 
Heisenberg model is recovered independently of p, and the system is hence ordered at T = 0. 
For A ^ 1 (but A 7^ 1) both limits p and p ^ 1 correspond to the 2D Heisenberg 
model with dilute bond-impurities. Thus the system should be ordered in these regimes as 
well. As p is increased from there will be an increasing fraction of singlets forming at 
isolated strong bonds in a background of weakly coupled spins. We argue that at a lower 
critical concentration p = pd, this leads to an order-disorder transition, in analogy with 
order-disorder transitions due to singlet formation in clean quantum antiferromagnets, such 
as the 2-layer Heisenberg model P| and various other dimerized models 0. As p is increased 
further, there must be another transition to an ordered state at p = Pc2, as the strong bonds 
start to dominate and the weak bonds effectively become impurities in a background of 
strongly coupled spins. As A is lowered the tendency to singlet formation diminishes, and 
one would expect the range [pci(A),]9c2(A)] to become smaller and eventually vanish at 
some A = A^m [0- Below we present numerical results supporting this picture, which is 
illustrated by the phase diagram outlined in Fig. 1. The solid circles at the phase boundary 
are results of our quantum Monte Carlo simulations, which are discussed below. 

We have used a modification of Handscomb's quantum Monte Carlo technique |ll],|12 



and averaged over 50-300 realizations of the random couplings in order to obtain results 
useful for extrapolation to the thermodynamic limit. We have studied systems of spins 
with periodic boundary conditions. In order to obtain ground state results for L = 4 — 10, we 
have carried out simulations at inverse temperatures (3 = J /T as large as 128, which for these 
system sizes is enough for all calculated quantities to have saturated to their T = value. 
A theorem by Lieb and Mattis []T^ guarantees that the ground state of a finite system with 



an even number of spins is a singlet, as long as all couplings are antiferromagnetic. We have 



therefore restricted the simulations to the subspace with zero magnetization (J^iS- = 0). 
We have also studied lattices with L = 32 at higher temperatures. In these simulations 
Monte Carlo moves changing the total magnetization were included. 

The sublattice magnetization m for a finite system can be defined according to 

= 35(7r,7r)/L2, (2) 

where S{n, n) is the staggered structure factor 

5(7r,vr) = -l^e-"(^V^^.)(^^^^^). (3) 

For the clean 2D Heisenberg model, spin-wave theory gives the leading size dependence of 
as ITSl 



m\L) = m\oo) + kL'\ (4) 

In Fig. 2 we show results for at a strong-bond concentration p = 1/4 for system sizes 
L = 4, 6, 8, 10, and various values of the disorder strength A. An approximately linear 
dependence on 1/L is seen for A < 0.8. We therefore fit straight lines to those points and 
extrapolate to obtain the sublattice magnetization for the infinite systems. For A = we 



obtain m = 0.276 ± 0.004, which is slightly lower [|1J] than the spin-wave result m = 0.303 
| 15| . For A = 0.70 the extrapolated is zero within statistical errors, and this should 



therefore be close to the critical disorder strength for p = 1/4. For A > Ac the scaling with 
system size must change to (l/L)"^ for large L, as ^(Tr, tt) saturates. We have also performed 
simulations at p = 1/8 and p = 1/2, and obtained Ac ~ 0.75 and Ac ~ 0.80, respectively. 
These points are shown as the solid circles in Fig. 1. The rest of the phase boundary is 
outlined schematically. The result for Ac(p =1/2) indicates that the critical concentration 
Pc2 as A ^ 1 is larger than the percolation threshold [p!0[| . 

For a given p, with Pci < P < Pc2, as A ^ Ac the spatial correlation length diverges 
as 6^'^, where 6 = |A — Ac|. The correlation length in imaginary time, C,t, diverges as 6~^'^, 



where z is the dynamic exponent |Tg]. With A > Lorentz-invariance is broken, and one 



expects 2; 7^ 1. The dynamic exponent can be determined by comparing the size-dependence 
of the staggered structure factor S{n,7r) and the staggered susceptibihty x(7r,7r). For a 
zero-temperature quantum phase transition the exponent rj for the algebraic decay of the 
spatial correlation function C{r) is defined by [|T^ 

where D is the spatial dimensionality. The staggered structure factor therefore diverges as 
^-u{2-z~rj)_ The staggered susceptibility is given by 

1 ^ 

= -.Y.^'^-^''-'"^ j dT{s]iT)sm), (6) 

and diverges as S~'^^'^~''i\ 

Finite-size scaling theory [|T8| gives the size dependence of ^(Tr, vr) and x(7r, tt) at the 
critical point: 

^(7r,7r)~L2-^-^ (7a) 

xin,n)^L'~^. (7b) 

Hence, if S'(7r, tt) ~ L^s ^nd x(7r, tt) ~ L"'^, the dynamic exponent is given by z = 'yx~ls- Iii 
Fig. 3, ^[^(vr, vr)] and ln[x(7r, tt)] are graphed versus ln(L) for two points which we estimate 
to be close to the phase boundary in the {p, A)-plane. Least-squares fits of straight lines 
give the slopes 7, = 1.01 ± 0.01, 7^ = 2.88 ± 0.09 for p = 1/4, A = 0.7 and 7, = 1.00 ± 0.02, 
'-f^ = 3.12 ± 0.15 for p = 1/2, A = 0.8. Hence we have two independent estimates for the 
dynamic exponent; z = 1.87±0.10 and z = 2.12±0.15. The indicated errors are of course 
only statistical. There are also errors, which we believe smaller than the statistical ones, 
due to the points investigated not being exactly on the phase boundary. 

The value of z has consequences for the behavior of the uniform susceptibiliy Xu at the 
transition point. The hyperscaling theory developed by Fisher et al. gives fl^ 



Xu = ^Y.^StS^^)^5'^^^-\ (8) 



Hence, depending on 2;, the uniform susceptibility diverges, remains finite or vanishes at 
the critical point. For "dirty bosons". Fisher et al. argued that the total compressibility 
(which corresponds to Xu) remains finite, and z = D ||T^. The spin model considered here 
can be mapped onto a hard-core bose system with particle-hole and SU{2) symmetries. 
The additional symmetries might in principle change the universality class from the one of 
the systems considered in Ref. [|17|. Our results for the dynamic exponent are, however, 
consistent with z = 2 = D. The results displayed in Figs. 2 and 3 indicate that 77 is close 
to —1, which satisfies the bound 77 < jl^. In principle, one can use finite-size scaling for 
A > Ac to obtain the correlation length exponent u. We do not presently have sufficient 



data for a precise estimate, but the bound u > 2/ D = 1 [|T3,0] appears to be satisfied. 

We have calculated the uniform susceptibility for L = 32 systems at temperatures T/J = 
0.2 — 1.0. Results for p = 1/2 and various A are displayed in Fig. 4. The disorder clearly 
enhances the low-temperature susceptibility. In view of the fact that Xu is non-zero at 
T = for the clean 2D Heisenberg model, it is unlikely that the susceptibility of the random 
systems vanishes as T 0. According to Eq. (^, this implies that z > 2, consistent 
with the above estimates from finite-size scaling. If z = 2, Xu approaches a constant as 
T at A = Ac. In ID, random exchange leads to a low-temperature divergence of the 
uniform susceptibility [^]. This is also predicted in higher dimensions for systems with 



longer-range interactions ||22[. The natural interpretation of this behavior is that some of 
the spins are effectively isolated from the rest of the system due to their local environment 
of strongly coupled spins. One would expect this behavior in the disordered phase of the 
model considered here as well, but not on the phase boundary if indeed z = 2. 

A non-zero T = uniform susceptibility implies that the disordered phase is gapless. In 
order to further investigate the spectrum, we have calculated the imaginary-time correlation 
function C(r) = l/NJ2i{Si{T)Si{0)), and used the maximum entropy analytic continua- 
tion procedure |2y] to obtain the real-time wave-vector integrated dynamic structure factor 
S{uj) = l/NJ2qS{q,u!). We have calculated S{uj) for both clean and random systems. In 
addition, in order to test the method, we have studied the case where the strong bonds are 
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arranged in a regular staggered pattern such that every spin belongs to a pair connected by 
a strong bond (p = 1/4). In this case one expects a gap for A larger than a critical value. 
Fig. 5 shows low-temperature results for L = 10 systems. For the staggered and random 
cases A = 0.8, and the ground states of the corresponding infinite systems are disordered. 
The staggered system exhibits a clear gap; S{uj) is a narrow peak centered at cu/J ~ 1 + A, 
corresponding to the energy required for a singlet-triplet excitation of a spin pair connected 
by a strong bond. For the clean system there is a broad maximum around u/J ~ 2, and 
a narrow peak close to u = 0. In the thermodynamic limit, long wavelength fluctuations 
of the order parameter lead to a 5-function peak at = 0. This peak is here shifted to a 
non-zero energy by the small gap present in the finite system. For the random case the peak 
at C(j = is due to localized, gapless excitations. There is also more weight at low energies 
than for the clean system, which we associate with excitations involving primarily the weak 
bonds. 

It would clearly be interesting to study the properties of a random quantum antifer- 
romagnet with a gap. Two coupled layers, each described by the hamiltonain (|^), should 
have both gapped and gapless disordered phases, depending on p, A and the (non-random) 
inter-layer coupling. Work on this model is in progress. 
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FIGURES 



FIG. 1. Proposed phase diagram of the random exchange model (|1]) in the {p, A)-plane. 
The solid circles are Monte Carlo estimates of transition points between the antiferromag- 
netic (AF) and gapless disordered (GD) phases. The curve is a schematic outline of the rest 
of the phase boundary. 

FIG. 2. The sublattice magnetization squared versus the inverse system size at a strong- 
bond concentration p = 1/4 and A = (solid circles), A = 0.5 (open circles), A = 0.7 
(solid squares), A = 0.8 (open squares), and A = 0.9 (solid triangles). Where not shown, 
statistical errors are smaller than the size of the symbols. The dashed and solid lines are 
fits to the A = and A = 0.7 data, respectively. 

FIG. 3. Finite-size scaling of the staggered structure factor (open circles) and the stag- 
gered susceptibility (solid circles) for p = 1/2, A = 0.8 (top) and p = 1/4, A = 0.7 (bottom). 
The straight lines are least-squares fits to the points. 

FIG. 4. The uniform spin susceptibility versus temperature for systems of size 32 x 32, at 
a strong-bond concentration p = 1/2. Solid circles are for A = 0, open circles for A = 0.5, 
solid squares for A = 0.7, open squares for A = 0.8, and solid triangles for A = 0.9. 

FIG. 5. Wave-vector integrated dynamic structure factors for L = 10. The dashed curve 
is for a clean system at /? = 32, the solid curve for a random system with p = 1/4, A = 0.8 
at P = 64, and the dotted curve for a staggered strong-bond pattern with A = 0.8. 
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